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We derive semiclassical laser equations valid in all orders of nonlinearity. With the help of a 
diagrammatic representation, the perturbation series in powers of electric field can be resummed 
in terms of a certain class of diagrams. The resummation makes it possible to take into account a 
weak effect of population pulsations in a controlled way, while treating the nonlinearity exactly. The 
proposed laser theory reproduces the all-order nonlinear equations in the approximation of constant 
population inversion and the third-order equations with population-pulsation terms, as special cases. 
The theory can be applied to arbitrarily open and irregular lasers, such as random lasers. 
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I. INTRODUCTION 

Interest in random lasers with coherent feedback [1, 2] 
and lasers based on chaotic microresonators without mir- 
rors [3] revealed a number of shortcomings of conven- 
tional laser theory [4, 5] that complicated its application 
to such systems. Among the properties that characterize 
random and chaotic lasers are strong openness, irregular 
spatial dependence of the refractive index, and, possi- 
bly, nontrivial shape of the resonator. Lasing in these 
systems can be accompanied by strong coupling between 
the modes, which requires a more careful treatment of 
nonlinear effects than is necessary for regular lasers. 

An essential part of a laser description is the choice 
of an appropriate basis of normal modes in which elec- 
tromagnetic field and other system functions can be ex- 
panded. In an open system it is not possible to de- 
fine a Hermitian eigenvalue problem whose eigenfunc- 
tions would form an orthogonal basis. Instead, one has to 
introduce a biorthogonal system of so called quasimodes 
as left and right eigenfunctions of a non-Hermitian op- 
erator. A number of methods to construct quasimodes 
were discussed in the literature. Among the earliest are 
the Fox-Li modes [6-8] which are useful in resonators 
with a preferred propagation direction and clearly de- 
fined transverse plane. In a more general setting, there 
were attempts to use solutions of the wave equation sat- 
isfying outgoing boundary conditions at infinity (Siegert- 
Gamow boundary conditions), with complex eigenfre- 
quencies corresponding to scattering resonances [9-11]. 
However, these modes diverge at infinity, which makes it 
problematic to use them as a basis [9], while some ways 
around this problem have been discussed in Refs. [10, 11]. 
Another possibility is to use the so called system-and- 
bath type approaches [12, 13], where cavities are de- 
scribed by an orthogonal system of wavefunctions of a 
Hermitian problem with an independent sets of modes 
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introduced for outside of the resonator. The openness of 
the cavity in this approach is reproduced by introduc- 
ing coupling between discrete set of inside modes and 
the continuous spectrum of outside modes. This way, 
the modal expansion of the field becomes possible both 
inside and outside of the cavity. 

In the present work we use a different approach based 
on keeping the spectral parameter of the outside outgoing 
field real, while making inside field to satisfy continuity 
conditions at the boundary of the cavity. This approach 
also results in the discrete spectrum of the cavity field 
with complex-valued frequencies, but, since the outside 
field is forced to depend on real spectral parameter, it 
does not diverge and is characterized by a constant flux. 
This approach was used in Ref. [14] for special kind of vi- 
bration problems and was adapted for lasers in Ref. [15], 
where respective modes were dubbed the "constant-flux" 
(CF) modes. One can introduce two adjoint biorthogo- 
nal systems of CF modes, which can be used to represent 
field inside the cavity. 

In the standard semiclassical laser theory [4, 5] las- 
ing modes are usually taken to coincide with quasimodes 
of the respective cavities, while their amplitudes and fre- 
quencies are found from equations based on perturbation 
expansions containing terms linear and cubic in the field. 
In random lasers this picture needs to be revised. First, 
it was shown [16, 17] that normal modes in the presence 
of gain differ from the passive modes even in the lin- 
ear approximation if the refractive index and/or unsat- 
urated population inversion are nonuniform. Second, it 
was realized [18-20] that self- and cross-saturation coeffi- 
cients before the cubic terms can have different statistical 
properties in different systems, leading to different mode 
statistics. Finally, it was pointed out [15, 21] that non- 
linear effects can significantly contribute to modification 
of the lasing modes compared to those of the empty cav- 
ity. A theory suggested in Refs. [15, 21] allowed for self- 
consistent calculations of not only lasing frequencies but 
also of the spatial distributions of the respective modes. 
Neglecting pulsations of the population inversion the au- 
thors of Refs. [15, 21] were able to derive equations for 
field amplitudes and frequencies beyond the usual third- 
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order approximation. 

In the present work we also generalize the conventional 
laser theory, but, unlike the approach of Refs. [15, 21], we 
do not begin by introducing a special approximation for 
population inversion. Instead, we carry out the perturba- 
tion expansion up to the infinite order in the field keeping 
all the terms which do not have fast temporal oscilla- 
tions. This also includes a part of population-pulsation 
contributions which is consistent with the slowly- varying- 
envelope approximation. The classification of the result- 
ing terms becomes possible due to a special diagram tech- 
nique developed to represent the terms of the expansion. 
This technique, however, differs significantly from usual 
Feynman diagrams widely used in solid state and high- 
energy physics, because we have to deal with terms of 
ever increasing degree of nonlinearity. Thus, diagrams in 
our technique are not used to literally represent each term 
of the expansion, but play a more limited role, as a tool 
assisting in the classification of the terms. Nevertheless, 
this technique allows for standard division of diagrams 
into connected and disconnected, with disconnected dia- 
grams submitting to easy resummation in terms of only 
connected ones. The connected diagrams can be classi- 
fied according to the order of magnitude of the popula- 
tion pulsations. The resulting laser equations generalize 
the nonlinear equations of Refs. [15, 21] in two respects. 
First, equations derived in this paper are dynamic, al- 
lowing for studying time dependence of the amplitudes, 
while the equations of Refs. [15, 21] can only describe sta- 
tionary lasing output. Second, our equations incorporate 
terms responsible for the population-pulsation contribu- 
tion in all orders of the perturbation theory; the equa- 
tions of Refs. [15, 21] are reproduced in our approach if 
only the lowest-order diagram is taken into account. 

The structure of our paper is as follows. In Sec. II we 
recall the definition and properties of the constant-flux 
quasimodes of open system. Standard semiclassical laser 
equations are written in Sec. Ill in frequency representa- 
tion for later convenience. Coupled equations for electric 
field, polarization, and population inversion are reduced 
to equations for the field alone in Sec. IV using infinite- 
order perturbation theory. In Sec. V we formulate the 
diagrammatic technique and resum the perturbation ex- 
pansion in terms of connected diagrams. In Sec. VI we 
reproduce the results of linear theory, third-order theory 
with population-pulsation terms, and all-order nonlinear 
theory in the constant-inversion approximation. Finally, 
we write corrections to the all-order theory that arc of 
the first order in the population pulsations. 



an electric field E(r, t) is governed by the wave equation 



e(r) J^E + V x (V x E) 



0. 



(1) 



where Gaussian units with the velocity of light in vacuum 
c = 1 are used. 

In the absence of gain and absorption the field will 
decay in time due to the openness. It is convenient to 
represent the decaying field as a superposition of certain 
normal modes, the quasimodes, that have only outgoing 
components outside the system. These modes can be 
constructed as families of constant-flux (CF) states [15] 
i/> k (r,u>) depending on a real continuous parameter ui. 
Explicitly, the CF modes satisfy the differential equation 



Vx[Vx^)]= w 2 ^( u ) 



(2) 



in the exterior of the cavity with the outgoing-wave 
boundary conditions at infinity. Inside the system, the 
same state satisfies a different equation: 



=V x 



V x 



fi 2 fc ( w )^ fc H. (3) 



For each w, the complex eigenfrequency fife(^) is quan- 
tized, as the eigenfunctions are required to match 
smoothly at the interface. 

Conjugate wavefunctions (f> k (r,uj) obey Eq. (2) out- 
side with the incoming-wave boundary conditions and 
the equation 



1 



=V x 



V x 



[o 2 fe H]*0 fc H. (4) 



inside the system. The CF functions and their conju- 
gates are biorthogonal and can be chosen to satisfy the 
condition 



L 



dr<j>* k (r 7 uj) ■ ip k ,(r,u) = S kk > 



(5) 



where the integration is over the interior I. 

A Fourier component E w (r) of the internal field can be 
expanded in the CF modes as 



E w (r)=e- 1 /2 (r) ^a fe ( W )^ fe (r, W ), 

k 

a k (w) = J dr Vefr) </>;i(r,a;) • E w (r). 



(6) 
(7) 



II. CONSTANT-FLUX QUASIMODES 

We consider an open system defined by real dielec- 
tric constant e(r), with e — 1 outside of the system's 
boundary. In the Coulomb gauge V • [e(r) E(r, i)] = 0, 



When continued to the exterior, this expansion yields a 
wave at the frequency u> propagating in the free space 
away from the system. In a stationary lasing regime u> 
becomes subjected to an equation, which has a discrete 
set of solutions corresponding to the frequencies of lasing 
modes. 
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III. SEMICLASSICAL LASER EQUATIONS 



where the pump Ano(r) is assumed constant in time. 



In the semiclassical theory of lasers [4, 5] the fields are 
described classically at the level of Maxwell equations 
and the active medium is treated by quantum mechan- 
ics. To this end, the wave equation (1) is written with a 
source term, the polarization P(r, t) of the gain medium: 

e(r) |e + Vx(VxE) = ~ 4 ^|i P ( r ' *)■ ^ 

In the simplest model, the active medium is a collection of 
homogeneously broadened two-level atoms. Their state 
is fully described by P(r,i) and the population-inversion 
density An(r,i) (difference between populations of the 
upper and lower levels per unit volume) . These functions 
satisfy the equations of motion [5] 

W + 2l± §t +l/2 ) P = - 2 ^ E ( r '*) A «M), ( 9 ) 



d_ 

dt 



d 



An-7||[Ano(r,t) - An] = jr;E(r,i) • ^P(M) 



dt 



(10) 



where d is the magnitude of the atomic dipole matrix 
element, v is the atomic transition frequency, and 7j_ 
(7|| ) is the polarization (population-inversion) decay rate. 
If the right-hand side of Eq. (10) vanishes, An relaxes 
to the unsaturated population inversion Ano(r, t), which 
is a measure of the pump strength. The coupled equa- 
tions (8)-(10) determine, in principle, electric field in the 
system, if Ano(r, t) is given. 

It is convenient, at this stage, to rewrite the equations 
of motion in the frequency representation. The Fourier- 
transformed variables are given in the form 



(11) 
(12) 

(13) 



An(r,i) 



-Re 




L 


7T 




= -Re 








lo 






~ 2tt , 


1 —00 



dwAn aJ (r) e 



that facilitates application of the rotating-wave approx- 
imation. Additionally, we assume that E w and P w are 
negligible outside of a small vicinity of v. Then the fre- 
quency squared can be approximated as u 2 = [y + uj — 



2v(uj — v) and the lasing equations become 



effectively the first-order differential equations in time, 
- e(r) {-v 2 + 2vuj) E u |Vx(VxE u )= 4tw 2 P w , 

(14) 

d 2 f°° 

(15) 

{-iui + 7||)An w = 27T7|| An (r) 6(u) 
i f°° 



IV. ALL-ORDER PERTURBATION THEORY 



A. Expansions for polarization and population 
inversion 



Equations (14)-(16) can be reduced to an equation for 
the electric field alone using perturbation theory in the 
field amplitude. In particular, one needs to construct an 
expansion of P u (r) in the (odd) powers of the field using 
Eqs. (15) and (16). Then this expansion is substituted in 
Eq. (14) producing the required equation for the field. In 
the conventional laser theory [4, 5] P w (r) is expanded up 
to the third order in E w (r), which yields the saturation 
terms in the rate equations. In this paper we will carry 
out the expansion up to an arbitrary order in the field's 
amplitude and use diagrammatic method to sort out the 
respective terms. 

We begin by neglecting the quadratic terms in Eq. (16) 
which gives the zero-order expression for the popula- 
tion inversion: Ani^ = 27rAno(r) S(u>). Substitut- 
ing this expression in the polarization Eq. (15), we 

obtain the first-order term for polarization Pj 
-i(d 2 /h-f ± )D(uj) An (r) 



,(i) 



I E w , where 



D{u>) = 



-1 -1 



l-i- 



7-L 



(17) 



This expression is substituted back in Eq. (16), from 

(2) 

which the second-order correction to the inversion An^ 
is determined. This iteration procedure yields the pertur- 
bation series for polarization and population inversion, 



P„(r)= J2 P ^(r), 

q odd 

An w (r)= Y, AntT^r), 

q odd 



(18) 
(19) 



in odd and even powers of the electric field, respectively. 
Henceforth we restrict the calculations to the case of 
scalar field. The general terms of these expansions are 
derived by induction in Appendix A. The resulting ex- 
pressions are 
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=2»fry„ ?rg _ 1 An (r)£>(o;) j dwi • • • du^E^E^ ■ ■ ■ E Uq _,EZ q _E„. 



W g -2+Wq-l 

x D|| (wi - cj 2 ) -D|| (wi - cj 2 + w 3 - w 4 ) • • • D\\ (loi - lu 2 H h w 9 _ 2 - 

x [-D(wi) + D*(w 2 )][-D(u;i - w 2 + w 3 ) + £>*(^2 - wi + w 4 )] • • • 

x [D(wi -0J2+0J3 + u q - 2 ) + D*(oj 2 - wi + cj 4 w 9 _ 4 + Wg_i)], (20) 



x D|| (wi - W2) (wi - ^2 + - w 4 ) • • • Dp (wi - cj 2 H h CL) g _ 2 - w 9 _i) 

x [-D(wi) + D*(w 2 )][£)(a;i - w 2 + w 3 ) + £>*(u; 2 - wi + cj 4 )] • • • 

X [D(U>1 - UJ 2 + W 3 W 9 _ 3 + Wg_ 2 ) + D*(cl) 2 - Wl + LJ 4 CJg_ 4 + LOg-l)] 

X - W 2 + W 3 + bJ q ) + c.c.(u> -> -cj), (21) 



where g is odd. The notation c.c.(w — > -co) stands for the 
first part of the equation to which complex conjugation 
accompanied by change of the sign of oj was applied. We 
also introduced the following definitions 



A = - 



2fi 2 7j_7|| ' 



1-i- 



7|| 



(22) 
(23) 



The above expressions describe nonlinear (q > 3) correc- 
tions to polarization and population inversion, which are 
used below to obtain equations for the field amplitudes. 



B. Equations for electric field 



Equation for the amplitudes (7) of normal modes, 



2™ jf dre-V^r) #(r,w)P u (r), 



(24) 

follows from Eq. (14) and biorthogonality condition (5). 
The right-hand side is a sum of the infinite number of 
nonlinear corrections to the polarization and is, therefore, 
a functional of all amplitudes <Xfc(w). By definition, a las- 
ing mode is such a solution which, in the limit t — > 00, ap- 
proaches a purely harmonic form, oc exp(— iwkt). Here ujk 
is some frequency, which is determined self-consistently 
from Eq. (24). In the frequency domain these solutions 
are characterized by delta-functional frequency depen- 
dence oc 8(u> — cjfe). 

Generally speaking, lasing modes are different from the 
quasimodes of the passive system. Indeed, an attempt to 
search for the solutions of Eq. (24) in the form ak(u>) = 
afc 8(ui — uik) leads, after the frequency integration (20), 



to the equation 

[fifc(a>k) - Wfc] a k S(lo - uj k ) 
= J2 W ki 1 ,...,k q a k^*k 2 ---al q _ i a kq 

q odd ki,...,k q 

x S(uj - uj kl + uj k2 hWfc ! - w fe ), (25) 



where the coefficients Wj^l 



are functions of the fre- 



quencies u>k,^>kn ■ ■ ■ ,^k q - Clearly, most of the terms on 
the right-hand side contain the delta functions that are 
different from 8{uj — ujk) on the left-hand side. In the time 
domain, these terms would introduce oscillations with 
frequencies different from u>k- This means that, in gen- 
eral, no stationary lasing solutions exist unless, for some 
reasons, the terms oscillating at the "beat" frequencies 
are small and can be neglected. Usually the selection of 
the slowly changing contributions to Eq. (25) is done by 
leaving only terms with pairwise coinciding indices fcj, 
so that the respective frequency differences cancel out 
(more on this procedure can be found below). However, 
since the number of u>ki in this equation is odd, one will 
always remain with the expression S(cu — LJk'), where k' is 
the index of one of the remaining uncanceled frequencies. 

This problem can only be resolved by requiring that, 
in a given lasing mode, I, each contribution aik of the 
quasimode k oscillates at the same frequency, and the 
general solution for ak takes the form of 



(26) 



Amplitudes aik can be shown to obey the equation 



- cji] aik 



aik' 



(27) 



where the coefficients vjj£, depend on all frequencies and 

amplitudes. An explicit expression for Vull is given be- 
low. One can see from this equation that the lasing 
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modes are those combinations of the quasimodes that 
diagonalize the matrix V kk , , while lasing frequencies are 
real eigenvalues of this matrix [16, 17]. 

One has to realize, though, that Eq. (27) is not an or- 
dinary eigenvalue problem, since the matrix V kk , itself 
depends on the amplitudes aik- Unlike linear eigenvalue 
problems, that determine frequencies for which nonzero 
solutions for the amplitudes can exist, solving Eq. (27) 
one must be able to find the frequencies, as well as 
the field amplitudes. This is possible, because the re- 
quirement that the respective eigenfrequencies must be 
real provides an additional constraint on solutions of 
Eq. (27) [21, 22]. 

If matrix V kk , is calculated in the constant-inversion 
approximation (see below), Eq. (27) reproduces the main 
result of Ref. [23]. However, while the derivation of this 
equation in Ref. [23] is only valid in the strictly station- 
ary limit, the arguments presented here can be extended 
to a nonstationary case. Indeed, we can repeat the above 
arguments for a weakly nonstationary situation, requir- 
ing that the amplitudes aik be slowly changing functions 
of time. Formally, we replace the mode expansion (26) 
with afc(w) = Yli a ik{u — ui), where ai k (uj — toi) is as- 
sumed to be sharply peaked at u) = u>i. In this case we 
can transform Eq. (24) to the time domain by expanding 
Qk(u) as Cl k « f^(cj;) + {u> - u)i)Q! k {wi), where VL' k (uj) is 
the derivative of Vt k with respect to the spectral param- 
eter u). It was found in Ref. [24] that, at least in the case 
of modes of a disk resonator, this derivative is not small 
and must be taken into account. In nonlinear terms we 
simply replace aik(u — oji) — > irai k {t) 8{lo — lji), which 



amounts to neglect of time derivatives of the nonlinear 
corrections. The resulting equation is 

j-i [1 - fl' k (uji)]^ + [fifc(wj) - aik(t) 

k' 

which, in the stationary limit, coincides with Eq. (27). 
Note that V kk , is now a slowly varying function of time 
via its dependence on the amplitudes ai»k"(t). The cor- 
rection due to Q' k is a new term, which has not been dis- 
cussed in any of the previous treatments of lasing dynam- 
ics. While it does not affect the steady-state solutions, it 
might change their stability, and is, therefore, important 
for strongly open cavities. More detailed study of its role 
is outside of the scope of this paper and will be presented 
elsewhere. 

The polarization matrix V kk , in the slowly varying am- 
plitude approximation can be presented as 

V k {l k ](t)=2Trv J dve-\r) <f k (r,un) fa (r,w,) m (r, t), 

. ( 29) 

where we introduced the nonlinear susceptibility 
r]i(v,t) = Pi(r,t) / Ei(r,t) defined as the ratio of the 
slowly varying polarization amplitude Pi(r, t) and the 
field Ei(r,i) in the mode I. The expression for the sus- 
ceptibility is found from perturbation expansion for po- 
larization P(r,t), Eqs. (18) and (20), and is given by 



» ? ,(r,t)=2ift7||D(a; I )Ano(r) ]T A 3 ^ ^ \E h (r,t)\ 2 \E h (r,t)\ 2 ■ ■ ■ (r,t)| 2 

q odd li,...,l q 

x D\\(w h -ui 2 )D\\{u h -u h +uj h - u u ) ■ ■ ■ D\\{u h - uj h + ■ ■ ■ + u lq _ 2 - u lq _ t ) 
x \D{uj h ) + D*{uj l2 )][D{u) h -w h +uj h ) + D*(uj h - u h +w u )]--- 

x [D(uj h -uj h +uj h uj [q _ 3 + uji q _, 2 ) + D* (uj h -u h +uj u uji q _ i (30) 



Here the order of nonlinearity q, introduced in Eq. (18), 
determines the number of different indices k, which take 
values from 1 to N m , where N m is the number of lasing 
modes. The superscript "r" at the sum symbol specifies 
that the possible values of the indices are restricted by 
the resonance condition 

Uh -ui 2 +oj h u lq _ 1 + oj iq - uji = (31) 

which ensures cancelation of fast oscillating terms. In the 
absence of accidental degeneracies, this condition implies 
that each of the indices h,h, ■ ■ ■ ,l q must be equal to one 
of the indices I2, U, ■ ■ ■ , l q -i, I- This leads, in particular, 
to the appearance of absolute squares of the field in the 
first line of Eq. (30). Moreover, the index l q effectively 



drops out of the equation, since the amplitude Ei q must 
be equal to some other amplitude Ei t . It is assumed 
that the slowly varying field amplitudes arc expressed in 
terms of quasimode components ai k (t) of the respective 
Zth lasing mode using 

E t (r,t) =e- 1 / 2 ( r )^a ife (i)Vfc(r,wO- (32) 

k 

Substituting Eqs. (30) and (29) in Eq. (28) one obtains 
a closed system of dynamic equations for ai k valid to all 
orders in the field amplitude. 

One of the fundamental difficulties of the theory of 
lasers is that the number N m of lasing modes is a priori 
unknown and depends on the strength of the pumping 



and the spatial distribution of the electric field in the 
cavity. In Ref. [23] this value is determined by the num- 
ber of possible solutions of Eq. (27) with real frequencies 
at a given pumping strength. This approach, however, 
does not take into account stability of the found solu- 
tions, which can only be determined by considering the 
time-dependent Eq. (28). Using this equation one could 
start by assuming that N m is equal to the size of the 
basis of quasimodes, N^, and study their time evolution. 
Those Ei which do not correspond to real lasing solutions 
at given pumping would decay to zero, and the number 
of lasing modes would be determined a posteriori without 
the need for a prior knowledge of N m . This approach is 
not free of difficulties cither, because of possible multi- 
stable behavior and hysteresis. Analysis of these issues, 
however, is beyond the scope of this paper. 



V. DIAGRAMMATIC TECHNIQUE 

A. Diagrammatic representation of the 
perturbation series 

In this subsection we present a diagrammatic technique 
developed to classify different nonlinear terms in Eq. (30). 
It should be noted, however, that our diagrams, unlike 
diagrams of the field or many-particle theory, do not pro- 
vide one-to-one correspondence between different terms 
of Eq. (30) and elements of the diagrams. The role of 
the diagrams here is more limited: we use them to clas- 
sify different pairing possibilities for the lasing mode in- 
dices l\, I2, ■ ■ ■ , l q , I in the perturbation series (30). Nev- 
ertheless, as it is shown below, this technique allows for 
classification and partial summation of the classes of the 
terms in a manner very similar to traditional diagram- 
matic methods. Unlike pairing of vertices in traditional 
diagrammatic techniques, which reflects Vick's theorem 
for creation-annihilation operators or Gaussian statistics 
of respective random processes, the pairing procedure in 
the situation under consideration hinges upon the con- 
dition expressed by Eq. (31). The resonance condition 
guarantees the absence of the fast oscillating terms, and, 
hence, the validity of the slowly changing amplitude ap- 
proximation. 

To construct a diagram X^- of order q = 1,3,..., we 
place q + 1 vertices in two columns as shown in Fig. 1. 
The left vertices are labeled h, I3, . . . , l q and the right 
vertices are labeled Z 2 , 14, . . . , l q -i, I. The vertex I is dif- 
ferent from the other vertices, because there is no sum- 
mation over the index I in Eq. (30). After that, each 
vertex on the left is connected with exactly one vertex 
on the right. The index j = l,...,N q labels all distinct 
connection possibilities in an arbitrary order. To obtain 
all diagrams of order q, we can, first, connect the vertices 
by (q + l)/2 horizontal links and then reshuffle the ver- 
tices, say, on the left without cutting the links. Thus, the 
number of possible diagrams of order q is the number of 
permutations N q = [(q + l)/2]!. Diagrams for q — 1, 3, 5 



I 



q-2' 

V 



/ 



q- 1 



FIG. 1: Labeling of vertices in a diagram of order q — 1,3,.. 



(a) 



U* / 



X 



11 



(b) 
h 



(c) 



X 



31 



X 



32 



FIG. 2: First-order diagram (a) and third-order dia- 
grams (b, c). 



are shown in Figs. 2 and 3. 

Each diagram specifies a particular contribution to the 
series (30). The latter will be written in the form 



r]i(r,t) = 2ih-y\\kno{v) D{uJi)X h 



(33) 
(34) 



q odd 



where each represents a partial sum in (•••)' 
Eq. (30), in which pairs of indices are chosen to be equal 
to each other according to the links connecting respec- 
tive vertices in the diagram. For example, the first three 






A 54 



X 



55 



FIG. 3: Fifth-order diagrams. The dash-dotted lines are 
the cuts that split disconnected diagrams into connected dia- 
grams. 
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diagrams correspond to the following expressions: 



(35) 



*3°i =E l^i 2 (r,t)| 2 I>||(wj - ui a )[D{ui) + D*(ui a )], 

(36) 

Xl 2 =^\E h (^t)\ 2 2Re[D(uj l2 )]. (37) 



The restriction l 2 ^ I in the diagram X 31 excludes the 
term with l\ = l 2 = h = I, which enters X® 2 . In general, 
the terms with more than two indices equal belong to 
the diagram in which the links connecting these indices 
do not cross each other. Another example are the fifth- 
order terms with l 2 = h = I4 = h 7^ I, which enter X® 2 , 
but not X^. Expression for arbitrary X^- is given in 
Appendix B. 



B. Resummation of the diagrams 

A diagram is called connected if it cannot be cut by a 
horizontal line without cutting a link. For instance, the 
diagrams X^ lt X^, X® 2 , and X® 3 are connected, while 
the diagrams X° 2 , X® 4 , X® 5 , and X® 6 are disconnected. 
To simplify the notation, we ordered all connected dia- 
grams before the disconnected diagrams for given q. We 
will label connected diagrams as X° ■ with 



yO _ yO 

qj ~ iv 



l,...,N q 



(38) 



where N q (< N q ) is the number of connected diagrams. 
The horizontal cuts separate disconnected diagrams into 
one connected diagram containing the vertex I (denoted 
by an unfilled dot in the graphic representation) and sev- 
eral connected subdiagrams without such vertex. The 
latter subdiagrams will be denoted as X q j, j = l,...,N q , 
where q + 1 is the number of vertices in the subdiagram. 
In place of the vertex I these diagrams have a vertex with 
the index Z g+1 which runs over all lasing modes, as the 

other indices lj. For example, the diagram X® 5 consists 
of X^j and A 31 and the diagram X® 6 consists of X^ and 
two diagrams An, where 



^ii-^|^ 2 (r,i)| 2 2Re[i?( W;2 )], 
X 3 i= J2 \Ei 2 (r,t)\ 2 \E u (r,t)f 



(39) 



h # u 

x£>||(w, 4 -o;, a )[£>(a;«J + D*(a; Ia )] 2 . (40) 

A general expression for X q j is given in Appendix B. 
Note that X q j is of the order q + 1 in the electric field. 
Connected diagrams contain (q — l)/2 nontrivial factors 

Dn + 1. 



Our diagrammatic technique possesses the basic prop- 
erty that disconnected diagrams are given by products of 
their connected parts. Thus, we can write for our exam- 
ples 



A 55 — A ll A 31> 
^56 = X U (^ll) 2 - 



(41) 
(42) 



The multiplicativity is due to the fact that the resonance 
condition of the type (31) is fulfilled for each connected 
subdiagram (see also Appendix B). The multiplicativity 
property allows to express the series (34) in terms of the 
connected diagrams as 

/ N q \ 00 / N q 

*= EE^ E EE^ 

\ q odd j=l j m=0 \ q odd j=l 



2-~*q odd q^ 
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1 J2q odd Sj=l Xqj 

This resummation formula is the main result of our pa- 
per. 

VI. LIMITING CASES AND DISCUSSION 

To make the meaning of Eq. (43) more transparent 
we apply it in several well-known special cases. We 
start with the linear approximation in Sec. VIA and 
consider the effect of gain-induced coupling of passive 
modes [16, 17]. In Sec. VI B we reproduce semiclassical 
equations of the standard third-order laser theory [4, 5]. 
In Sec. VI C we discuss all-order nonlinear theory in the 
approximation of constant population inversion [15, 21 
23] and derive, using our theory, the first diagrammatic 
correction to it. 



(43) 



A. Linear gain-induced mode coupling 

In the linear approximation to the polarization Pi = 
rjiEi [Eq. (30)] only the lowest diagram X^ contributes 
to Xi. Equation (28), where the matrix Vuu, is calcu- 
lated using Eqs. (29), (33), and (35), yields the following 
equations for the slowly-varying amplitudes: 

E |^ fefe ' ^ + i \P kk ' - Skk'Ul] j aik> = 0, (44) 

where the fl' k term is henceforth neglected. The fre- 
quency matrix 

fifcfc'(w) = tt k {u>)6 kk > - Vfcfc'M, (45) 

is modified by the linear gain term, 

v kk ' M = 

- 2inv-^-D(u) £ rfr^M #(r, w) Mr, w), (46) 
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proportional to the overlap integral. Clearly, the ma- 
trices Vkk> and flkk' are nondiagonal if the pump or di- 
electric constant are not homogeneous. In this case the 
biorthogonal quasimodes of the system ip k and <f> k are no 
longer ip k and <p k , but are determined by the right and 
left eigenvectors of flkk' according to 

Mr,u J ) = J24k'^(r,u J ), (47) 

k> 

0;(r,w) = ^o^(r,w). (48) 
k> 

The right eigenvectors a\ T k ] are normal modes of Eq. (44) 
whose amplitudes ai (t) obey the equation 

-hi +i [TIi(uji) - bJi] ai = 0, (49) 

where f^(w) are eigenvalues of Qkk'i 1 ^ 1 )- We recall 
that electric field in the mode I has a time dependence 
ai(t) exp(— iujit). Thus, the lasing frequency ooi is deter- 
mined from the requirement 

Repi^t)] =lji (50) 

and the threshold condition for this mode is 

lm|n,(wi)]=0. (51) 

As follows from Eq. (49), the mode amplitudes diverge 
exponentially above the threshold. Hence, applicabil- 
ity of the linear approximation is limited to the pump 



strength below or at the threshold. However, the basis of 
normal modes can be used as a starting point in nonlinear 
theories. 



B. Third-order theory 

To obtain an approximation to Xi of the third-order 
in the field, we keep the diagrams X^ and in the 
numerator and the diagram An in the denominator of 
Eq. (43) and expand the latter: 

X^X^+X^X^+X^. (52) 

It is convenient to write lasing equations in the basis of 
quasimodes ip k and <f> k that diagonalize the linear part. 
Due to nonlinear effects, the lasing modes above the 
threshold, 

£7,(r,t) = e-V 2 (r)^a, fe (i)^(r,u; ; ), (53) 
k 

are, in general, linear combinations of individual quasi- 
modes. Equation for the amplitudes ai k (t) follows from 
Eq. (28), after taking into account the results of linear 
theory (Sec. VIA), and has the form: 



aik +i [fife(wi) - ui] a ik = (4— \ D(ui) f dr An ^^ - (j)* k (r,uji) Ei(r,t) 

«7|| \ n 1±J J i V e ( r ) 



n„ 



J2 \Ev (r, t)\ 2 {2Re[D(w,/)] + (1 - 6 lv )D\\ (w, - w v ) [D(ui) + D*(uj v )} } , k = 1, . . . , A b , (54) 
v=i 

I 



where A b is the size of the basis of quasimodes and 
A m is the number of lasing modes. The lasing frequen- 
cies loi and the mode thresholds need to be determined 
from these A b x A m equations in the stationary regime 
ai k = using, e.g., a selfconsistent procedure described 
in Ref. [21]. 



In some cases the standard assumption of a traditional 
lasing theory that the lasing modes coincide with the 
quasimodes of the cavity remain valid. In this case the 
total number of equations (54) is reduced to A b since the 
amplitudes are approximated as aik(t) = ai(t)Si k . Repre- 
senting a/(t) = \/Ti exp(it^) and separating the real and 
imaginary parts in Eq. (54) , we obtain 2 Ab real equations 



for the intensities /; and phases tp { , 



i l -21mp l (oj l )]I l 

2nv / d 2 y 

Tp l + Re[Q;(u;;)] -w; 
-kv ( d : 



I[Re 



D(uj0J2BwIi'{---} 



Im 



D(lui)J2Bw !,{■■■} 



, (55) 



(56) 



The terms enclosed in the braces are the same as those 
in Eq. (54). We defined overlap integrals for the quasi- 



9 



modes as 



Bw = 



where 



dr 



An (r) _* 



{e(rW 



4>i (r,wj)^/(r,a;j) 



(57) 



The lasing frequencies u>i are determined, together with 
the stationary intensities, from Eqs. (55) and (56) in the 
stationary regime = and lp l = 0. The rate equa- 
tions (55) and frequency equations (56) generalize the 
standard third-order semiclassical theory with the self- 
and cross-saturation terms [4, 5] to systems with strong 
openness and arbitrary distribution of refractive index. 



C. Constant-inversion approximation and 
corrections 

In many physical situations the population inversion 
An(r,t) is approximately constant in time. Oscillations 
of the population (population pulsations) are responsible 
for the terms proportional to D\\{wi — lui>), I ^ in 
the expansion of the nonlinear susceptibility (30) and, 
hence, in the expressions for the diagrams. Comparing 
the terms in the braces in Eq. (54), we can conclude that 
the population pulsations can be neglected if D»(cji — 
cui>) -C 1, i.e., \uji — ^> 7|| for all lasing frequencies 
UJi ^ u>i>. 

In the approximation of constant inversion only the 
diagrams X^ and In, which do not contain the D\\ 
functions, contribute to Xi (43). Equation (28) for the 
coefficients of expansion (32) of the electric field becomes 



aik + i[Vt k (u){) - Lo{\aik = 2ttv 



d 2 



dr 



An (r) 



1 + 27F 7 _l- 



Y v \Ev{v,t)\nMD^v)\ 



(58) 



In contrast to Eq. (54), the linear contribution here is 
not diagonalized and is contained in the right-hand side. 

Equation (58) is valid in all orders in nonlinearity if 
the population inversion is constant. With the help of 
Eq. (43) it is straightforward to write out the correc- 
tions due to the population pulsations. The terms of the 
first order in Dn <C 1 are contained in the diagrams 
and X 3 i [Eqs. (36) and (40)], so that X\ can be approxi- 
mated as 



Xi 



x° 1 + x? 



31 



1-Xu -X: 



31 



(59) 



These corrections modify both the numerator and de- 
nominator of Eq. (58), which can now be presented as 



d 2 f 

fryi. Jx 



An (r) 



#(r,w,)£iO-.t)(l-Tn) 



1 + 2l4^ ^1' l^'( r ' t)\ 2 2Be[D(ui>)] (1 - T d ) 



(60) 



Tn=2^El^(r,*)l^||(«»-««0 

x [£)(w,) + D*(wi/)], 

Td=2^El^(r,i)| 2 J D||(^-^) 

[DM + D*^)] 2 
2Re[D{u l ,)] 



(61) 



(62) 



Taking into account that the electric field in the lasing 
modes has a typical magnitude of H^j±j^/d, one can see 
that deviations from the constant-inversion approxima- 
tion is of the order of D\\{ui — ujv) [the term D{loi) is of 
the order of unity since frequencies of lasing modes are 
concentrated within the width of the gain profile]. The 
difference between lasing frequencies can be estimated as 
\oji — oji'\ w r y±/N nl . Then the condition Du <C 1 can be 
expressed as N m j^/^f± <C 1. Given that 7|| is usually sev- 
eral orders of magnitude smaller than j±, this condition 
is in most situations fulfilled. It was reported in Rcf. [22], 
however, that nonlinear interaction between lasing modes 
can push their frequencies toward each other making the 
intermode spectral interval much smaller than the typical 
value given above. Such pairs of modes can result in sig- 
nificant corrections to the constant-population approxi- 
mation. Adding to the expansion additional connected 
diagrams with up to q + 1 vertices, one can improve the 
constant-population approximation by constructing las- 
ing equations valid in the order (q — l)/2 in Du. 



VII. CONCLUSIONS 

We presented a diagrammatic semiclassical laser the- 
ory valid in all orders of electric field. The original pertur- 
bation series in the powers of the field can be resummcd 
in terms of a certain class of diagrams, the connected 
diagrams. The resummation allows to construct a con- 
trolled expansion in the small parameter 711/7.1, which 
is a measure of population pulsations, while treating 
the nonlinearity exactly. Our lasing equations gener- 
alize the all-order nonlinear equations in the constant- 
inversion approximation and the third-order equations 
with population-pulsation terms. The use of constant- 
flux quasimodes as basis functions makes it possible to 
apply the theory to strongly open and irregular systems, 
such as random lasers and lasers with chaotic resonators. 



Acknowledgments 

Financial support was provided by the Deutsche 
Forschungsgemeinschaft via the SFB/TR12 and 
FOR557 (O.Z.) and by PSC-CUNY via grants No. 62680- 
00 39 and No. 61788-00 39 (L.D.). 



10 



Appendix A: Derivation of Eqs. (20) and (21) 

For a few lowest values of q the validity of Eqs. (20) 
and (21) can be checked directly. To prove these re- 
lations by induction, we assume that An£? -1 ^ is given 
by Eq. (21), then derive P^ q) (20), and, finally, ob- 
tain An£ 9+1) . 

According to Eq. (15), 



(Al) 



we first exchange the labels ui\ -H- u>2, ^3 W4, . . . , 
u;,j_2 — > Wg_i and then define the variable 

LO q -2 = UJ - Ld' - OJl + Ld 2 h Wq_3 + Ul q -\. (A4) 

Transforming the integral over to' as 

y dJ D|| (w - w') £' u ,'£' w _ w '_ Wl+W2 _... +W5 _3 +w<; _ 1 
= J duj q - 2 D\\(ui-u} 2 -\ hcj 9 _ 2 - 

X -E'W — Wl+OJ 2 ",-2+^,-1^,-2 (A5) 



Substituting An£ 9 ^ from Eq. (21), we immediately see 
that the factor before the integral in Eq. (20) is re- 
produced. To calculate the integral, we consider sep- we obtain the part of the integrand proportional to 

arately the two contributions: An^_^} = AnJ_J, + ^ ^ 2 Wl 

Next we derive Ani? using Eq. (16), 



An,.,,,., 



and 



An 



~ (9-1) 

, where An w is given explicitly by Eq. (21) 

(9-1)1 

is c.c.(cj — > — w). When integrating 



An,^ 1 ' 



D|,(a;) du' E* U ,_ U P$ 



(-/) 



7rfi7|i 

+ c.c.(w -> -w). (A6) 

Clearly, the factor before the integral in Eq. (21) follows 

uj q -i = to' — to + L0\ — uo 2 + • • • — w g _3 + u) q -2- (A2) after substitution of pffl (20). Introducing the new vari- 
able 



~ (9-1) 

An ul _ UJ , we introduce a new variable 



Then the integral over ui' becomes 

J dio 1 D,| (w - a/) S w ^*,_ w+u , i _ u , 2+ ..._ w<; _ 3+W(; _ 2 
= y dw 9 _iD|| (wi - u; 2 H h w 9 _2 - Wg-l) 

X W,_2+W,-l-E' 1 j,- 1 - (A3) 

Comparison with Eq. (20) shows that An^_ w , contribu- 
tion yields the part of the P^ integrand proportional to 



LO q = J - LOl + 0J 2 CJ 9 _2 + Uq-1 (A7) 

we rewrite the lj' integral 

y duj' D{u) r ) E*, -u 1 +uj 2 w g _ 2 +",-i 

du q D((Jl - LJ 2 H U)q-1 + U) q ) 



D{lj\ — 0J2 - 



Uq-2). When integrating 



~(9-l) 



X ^-W+Wl-W 2 H ^ q -l+Ul q ^^q ■ 



This completes the proof of Eqs. (20) and (21). 



(A8) 



J 



Appendix B: General expressions for diagrams 

The diagrams X^- (q = 3, 5, . . .) are given by the following analytical expression: 

X° qj = ]T \E h (r, t) | 2 \E h (r, t)| 2 • • • \E lq _ x (r, t) | 2 

x°. 

x D\\(ui h -u l2 )D\\{u h -ui h +uj h - lj u ) ■ ■ ■ D\\(u) h -uj l2 H V uji q _ 2 -W; 5 _J 

x [D(w, 1 )+D*(w, 2 )][£)(w, 1 -w h +w i3 ) + J D*(w/ 2 - u h + lj u )]--- 

x [£>(u; h -w, 2 +uj h uj lq _ 3 +u)i q _ 2 ) +D*(uj h -oj h +uj u uji q _ i +uj lq _ 1 )}. 



(Bl) 



The symbol denotes a summation over the lasing-mode indices fe, h,..., l q -i — 1, . . . , Am according to the rules: 



A' 



(i) if the indices k and lj are connected in the diagram X®j, set k = lj and (ii) the terms with four or more indices 
l\, I2, ■ ■ ■ , l q , I equal are excluded unless the links connecting the affected vertices do not intersect. 

We denote by A°- connected diagrams and subdiagrams containing the vertex I. With the ordering of A° such that 
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the connected diagrams come before the disconnected diagrams for a given q, we can identify = X^ (j < N q ), 
where N q is the number of connected diagrams. The subdiagrams X q j without the special vertex have a variable 
index l q+ i in place of the fixed index I. The analytical expression for these diagrams has the form 

X gj = £ \E l2 (r, t)\ 2 \E h (r, t)\ 2 ■ ■ ■ \E lq+1 (r, t)\ 2 

x qj 

x D\\{u h - uj h ) D\\{u h -Lu h +u h - u u ) ■ ■ ■ D\\{u h - uj h + ■ ■ ■ + u lq _ 2 - 0Ji q _J 
x [D(Lu h ) + D*(LOi 2 )}[D(Lu h -co h + lui 3 ) + D* (uji 2 -LJ h +uji 4 )]--- 

x [D^i-l - uji 2 + uji 3 — Ldi q _ 1 + Lui q ) + D*(lui 2 - o;; 1 + w; 4 uji q _ 2 + u>i q+1 )]. (B2) 



The sum £ is defined analogously to the sum £ above. 



x a 



X° 



Note that the diagrams X^ and X^ are of the order q— 1 
in the electric field, while the subdiagrams X q j are of the 
order q + 1. 

To show the multiplicativity property, let us assume 
that a disconnected diagram X^- can be cut in two, pos- 
sibly disconnected subdiagrams, X q ,j, and X q »j", such 
that q = q' + q" + 1. We will label the vertices of X q iiyi 
ash,..., lq"+i and identify them with the first q" + 1 ver- 
tices of the diagram • . The last q' + 1 vertices of X^ , 
l q "+2, ■ ■ ■ ,l q ,l, are also the vertices of X { q \y. We need to 
show that, in the sum £, the two groups of indices can 



x° 



be split between the sums £ and £ , respectively. 



E 

x", ., 
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This would mean that the arguments of the functions D 
and D|| can contain only the indices belonging to one of 
the groups. This is, indeed, the case, since the arguments 
having less than q" + 1 frequencies contain only the in- 
dices from the first group, while in the arguments with 
more frequencies the first q" + 1 frequencies cancel due 
to the cut as 



L0 2 + 0J4 H 1- 0Jq"+l- 



(B3) 



According to this equality, the number of nontrivial fac- 
tors D|| ^ 1 in Eqs. (Bl) and (B2) is (q- l)/2 minus the 
number of cuts (in a disconnected diagram) . 
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